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ABSTRACT 

We carry out iV-body simulations to examine the effects of dynamical interactions on 
planetary systems in young open star clusters. We explore how the planetary pop¬ 
ulations in these star clusters evolve, and how this evolution depends on the initial 
amount of substructure, the virial ratio, the cluster mass and density, and the initial 
semi-major axis of the planetary systems. The fraction of planetary systems that re¬ 
mains intact as a cluster member, /bpSi is generally well-described by the functional 
form /bps = /o(l + [a/ao]°)~^, where (1 — /o) is the fraction of stars that escapes from 
the cluster, ag the critical semi-major axis for survival, and c a measure for the width 
of the transition region. The effect of the initial amount of substructure over time t 
can be quantified as /bps = A{t) + B{D), where A{t) decreases nearly linearly with 
time, and B{D) decreases when the clusters are initially more substructured. Provided 
that the orbital separation of planetary systems is smaller than the critical value oq, 
those in clusters with a higher initial stellar density (but identical mass) have a larger 
probability of escaping the cluster intact. These results help us to obtain a better 
understanding of the difference between the observed fractions of exoplanets-hosting 
stars in star clusters and in the Galactic field. It also allows us to make predictions 
about the free-floating planet population over time in different stellar environments. 

Key -words: Open clusters and associations: general - A-body simulations - extra¬ 
solar planets - planetary systems 


1 INTRODUCTION 

During the last decade, the Kepler space mission and other 
transit, radial velocity surveys, microlensing, and imaging 
surveys, resulted in the discovery of thousands of exoplanets 
and exoplanet candidates (e.g., Borucki et al.|20TT Batalha 


et al. 20131, demonstrating that exoplanets in the Galac¬ 


tic field are common. Combined with new developments in 
computational astrophysics, the prospects for deepening our 
understanding of the formation and dynamical evolution of 
exoplanets are promising. 

In contrast to the large number of exoplanets found in 
the Solar neighborhood, only a handfull of exoplanets have 


been detected in open clusters, such as Hyades ( Guenther et 
al.|2005 1, M37 (Hartman et al.|2009 1, NGC2158 (Mochejska 
at al.|2006|, NGC7 789 ( |Bramich fc Horne|2006|, NGC124 5 
I Burke et L.||2006 l and NGC7068 ( Rosvick fc Robb||2006 1, 
and in globular clusters, such as uu Centauri (Weldrake et 


al. 


20081 


47 Tucanae (Weldrake et al. 20051, NGC6378 
([Nascimbeni et al.||2012 and NGC6791 ( Mochejska at al. 
2005 Montalto at al. [20(37 1. Although the study of Meibom 
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et al. {20131 suggests that the planet frequency in open clus¬ 
ters is similar to that in the field, statistically stronger ob¬ 
servational evidence is necessary strengthen this argument. 

Observations have indicated that most stars in the 
Galactic field were formed in hierarchically structured star¬ 
forming regions ( Bressert et al.|2012 (, often (but not always) 
in embedded clusters (e.g., [Clarke et al.|2000[ |Lada fc Lada| 
20031 that may or may not be gravitationally bound (e.g., 
Kruijssen 20121. This may also be true for our own Sun, 


which is believed to have been born in a young star clus- 


ter with several thousands of stellar members (see 

Adams & 

Lauslilin|2001 

Portegies Zwart|20091 Adams|2010 

IDukes & 

Krumholz|2012 

Pfalzner 2013 Parker et al. 120141 for discus- 


sions). In such dense stellar environments, close encounters 
between planetary systems and the stellar populations sur¬ 
rounding them can result in orbital reconfiguration or decay 
of planetary systems (e.g., Spurzem et al. 2009 Boley et 


al.|2012 Hao et al.|2013 1 and in the subsequent production 


of free-floating planets. Although free-floating planets are 
expected to be common in young star clusters and in the 
Galactic field (e.g., Pacucci et al. pol^ , only few have been 
detected through microlensing campaigns (e.g., [Sumi et al!] 
2011 Di Stefano|2012 Gaudij2012 1, and deep imaging sur- 































































































2 Zheng, Kouwenhoven & Wang 


veys of young star clusters (e.g., Lucas et al.|2006 Quanz et 
al.||2010 Delorme et al.|2012 |. 

The currently known planetary systems and free- 
floating planets are unlikely to have formed at the places 
where they are found today, and the orbital properties of 
these planetary systems may have been affected by their 
internal dynamical evolution and by gravitational interac¬ 
tions with the environment in which they were born. To 
understand the observed orbital properties of exoplanets, it 
is necessary to consider the dynamical evolution these sys¬ 
tems have undergone since their formation, prior to, during, 
and after the dynamical process that led to the escape of 
these planetary systems from their natal environment. 

Computer simulations are crucial for developing a good 


Table 1. Initial conditions for our reference model. 


Parameter 


Default value 


Number of stars 
Number of planets 
Structure 
Virial state 
Stellar mass function 
Half-mass radius 
Planet mass 
Semi-major axis 
Eccentricity 
Orbital orientation 
Orbital phase 
Simulation time 


N„ = 1000 

Np = 1000 

D =1.6, 2.3, 3.0 
Q =0.3, 0.5, 0 .7 


Kroupa 120011; O.OSMq ^ M ^ IOMq 
1 pc 

Mp = 1 Mj 
a = 100 AU 
e = 0 
Random 
Random 
t = 50 Myr 


planetary systems. Many studies (e.g., Laughlin & Adams 


1998| Bonnell et al.||2001| Smith & Bonnell||2001| |Davies & 

Sigurdsson 2001 

[Adams et al. 2006| h'regeau et al.| 2006 the work of Parker & Ouanz 

(20121, and our initial con- 

Liu et al. 2013 

Hao et al. 20131 have focused on the dif- 

ditions resemble theirs. In this paper, however, we study 


ferent aspects of planetary dynamics in star clusters. In 


an extensive theoretical study, Spurzem et al. (20091 de¬ 


rive collisional cross-sections for changes in the orbital ele¬ 
ments during encounters between stars and star-planet sys¬ 
tems in star clusters. Previous studies employing direct N- 
body simulations of single-planet systems in clustered en¬ 
vironment have mostly focused on open and globular clus¬ 
ters that are initially in virial equilibrium and have smooth 
initial density profile (e.g.. Hurley fc Shara||2002 Spurzem 
|et al.|[^09[ ). Observations, however, indicate that star clus¬ 
ters are preferably born with a high degree of substructure 
(e.g., Cartwright fc Whitworth|2004 Schmeja|2011 1 and of¬ 


ten in a subvirial state (e.g., [Peretto et al.||2006[ |Proszkow 


et al. 20091. Inspired by the earlier works of Goodwin & 


Bastian (20061, Parker & Quanz (20121, Craig & Krumholz 
( |2013P aiid |Pa!cucci et al.| ( |2013[ ), we carry out A'^-body sim¬ 
ulations of clusters with different initial morphologies and 
initial virial states, and a population of planetary systems 
with a wide range of semi-major axes. 

This article presents a theoretical study on the evolu¬ 
tion of planetary populations in star clusters, with a focus on 
how the dynamical evolution of these planets is affected by 
the initial properties of the star cluster and the initial semi¬ 
major axes of their orbits, and the subsequent dynamical 
evolution of the liberated population of free-floating plan¬ 
ets. In Sectionj^we describe our methods and assumptions. 
In Section we present the results of our simulations, and 
finally, we summarise and discuss findings in Section]^ 


2 METHOD AND TERMINOLOGY 


2.1 Initial conditions 


We use the NB0DY6 package (Aarseth 20031 to carry out 


A-body simulations of evolving star clusters that initially 
contain a large population of (single-planet) planetary sys¬ 
tems. Our initial conditions represent the state of the system 
shortly after most of the gas is expelled. The initial condi¬ 
tions are generated using the open-source package MCluster 
( Kiipper et al.|2011 l, which we modified to suit our require- 
ments. Table [H provides a summary of the initial condi¬ 
tions for our reference models. Our study was inspired by 


much broader range in planetary orbital separations. Unlike 
Parker & Quanz ( 2012| , we include an external tidal field, 


we normalise the initial size of the clusters to their half-mass 
radii (see the discussion below), and we do not include pri¬ 
mordial stellar binaries. All results described in our study 
represent the averages of twenty realisations of each model, 
unless stated otherwise. 

Our reference models are open cluster-sized systems 
composed of Nsp = Ns + Np = 2000 individual bodies, 
among which Ns = 1000 stars and Np = 1000 planets, each 
orbiting a star. Stellar masses are drawn from the |Kroupa| 


(20011 initial mass function in the range 0.08 — 10 Mq. 


The lower limit corresponds to the hydrogen-burning limit 
(e.g., Karttunen et al.|2003 l. Our choice for the upper limit 
is roughly consistent with Weidner & Kroupa (20041 and 


Weidner et al. (20131, although our random sampling of 


the masses from the IMF results in upper mass limits of 
8.6 ± I.IMq for the individual clusters (and it should be 
noted that the relation between the stellar upper mass limit 
and the cluster mass is still under debate, see, e.g., |Cervino| 


et al. 2013a|b I. Each planet is assigned a mass of 1 Mj 
(~ 9.5 xTWNUq). For the purpose of our study, the plan¬ 
ets can be considered as test masses, as their effect on the 
dynamics of the stellar population is negligible. The ini¬ 
tial total mass of the clusters in our reference model is 
Mci ~ 600 Mq , a typical value for open clusters. 

Although the vast majority of know exoplanets orbit 
their host star at much smaller separations, wide planetary- 


mass objects have been observed (e.g.. 

Lafreniere et al. 2008 

Marois et al.||2008 

||2010 

[Lagrange et al.||2009[ [Kraus et al.| 

2014 

Naud et al. 

2014 

Bailey et al.||2014 

1, and their exis- 


tence is also supported by numerical studies involving forma¬ 
tion at wide orbits or dynamical interactions in multiple star 
systems and planetary systems that lead to outward scatter- 


ing (e.g., Portegies Zwart & McMillan 

2005 

Stamatellos & 

Whitworth[2008[|Malmberg et al.|2011[ 

Boss 

2011[|Vorobyov 


reference model are assigned a planet with semi-major axis 
a = 100 AU. Much larger values of a result in frequent dis¬ 
ruptions, while single-planet systems which much smaller 
values of a are mostly dynamically inert in the stellar en¬ 
vironments we model. For a similar study involving plan¬ 
ets with a = 5 AU and a = 30 AU we refer to the work 















































































































































Planetary systems in young star clusters 3 


of Parker & Quanz (20121. Note that among planets with 


identical a, those orbiting more massive stars have a larger 
collisional cross section (due to gravitational scattering), but 
also have a larger binding energy that is proportional to the 
host star’s mass. To study the relationship between the sta¬ 
bility of planetary systems and the initial semi-major axis, 
a, we also study systems with a wide range of semi-major 


axes (1 AU ^ a ^ 10 000 AU) in Section 3.3 We merely 


use these initial conditions for the semi-major axis to study 
the dynamical evolution of such systems; we do not imply 
that they formed at these semi-major axes. Each planet is 
assigned a random orbital phase in a circular (e = 0), ran¬ 
domly oriented orbit. 

Using MCluster we generate star clusters with with 
varying degrees of substructure, quantified through the frac¬ 


tal (substructure) parameter D (see Cartwright & Whit- 


|worth|2004| [Goodwin fc Whitworth|2004| [Allison et al.|2009[ 
for details). We adopt three values, D = 1.6 (highly clumpy), 
D = 2.3 (intermediate substructure), and D = 3.0 (homo¬ 
geneous on large scales). The procedure for generating these 
substructured star clusters in MCluster is fully described 
in the appendix of |Kupper et ah' (20111, and briefly sum¬ 
marised below. Initially, a single star is placed at the origin 
of a box. The box is split up into eight smaller boxes, and 
in some of these a new star is placed; the probability that 
a star is placed in a box is where D is the fractal 

parameter. The boxes containing stars are again subdivided 
into eight boxes, and the procedure is repeated. During each 
step, velocities are drawn from a normal distribution centred 
on the parent body, and subsequently rescaled, so that the 
group of stars in the eight boxes are in virial equilibrium. 
After generation of a large number of stars, a subset of Ns 
stars within the unit sphere is drawn from the cloud of stars. 
Finally, the cluster positions are scaled to the desired radius 
of the cluster, and the velocities are scaled to obtain the 
desired initial virial ratio. 

The virial ratio Q of a star cluster is defined as Q = 
—K/P, where K and P are the total kinetic energy and 
potential energy, respectively. We study clusters that are 
initially in virial equilibrium (Q = 0.5), clusters which start 
in cool collapse (Q = 0.3) and clusters that are initially 
expanding (Q — 0.7). Each cluster is scaled to an initial 
(three-dimensional) half-mass radius rhm = 1 pc, which is 
typical for young star clusters in this mass range (e.g.,|Lada] 


& Lada 20031. Consequently, the radius of the sphere en¬ 
closing all stars equals 2}^^ pc « 1.26 pc. The star clusters 
are assigned Solar orbits around the Milky Way centre. The 
tidal radius Rt of a star cluster at any time is then given by 


R, 




1/3 


; 6.65 


Mci 

1000 Mq 


1/3 


pc 


( 1 ) 


Table 2. Classification of the planets in our simulations. 



Orbiting a star 

Free-floating 

Cluster 

Bound planetary 

Bound free-floating 

member 

system (BPS) 

planet (BFP) 

Escaped 

Escaped planetary 

Escaped free-floating 

from cluster 

system (EPS) 

planet (EFP) 


however, that for our choice of initial conditions, the virial 
radius depends on D. A star cluster that is homogeneous on 
large scales D = 3 has typically a larger virial radius than a 
highly substructured cluster (D = 1.6). Consequently, star 
clusters with smaller D tend to evolve faster. Numerical tests 
indicate that for our initial conditions the relation between 
physical time /Myr and N-body time /nb can be expressed 
as 

tMyr ~ tNB(0.37 -|- 0.13D) , (2) 

such that our total integration time of 50 Myr corresponds to 
roughly 93 N-body units for D = 1.6 and about 66 N-body 
units for D = 3.0. It should be noted, however, that this 
relation is approximate due to stochastic variations resulting 
from the relatively small number of stars in our simulations. 


2.2 Dynamical status of the planets 

In our analysis of the evolution of the planetary population 
we distinguish between four dynamical categories, based on 
whether or not a planet is gravitationally bound to a star, 
and whether or not it is a member of the star cluster. A 
planet is considered bound to a star when (i) the star and a 
planet are each other’s mutual nearest neighbours, and (ii) 
the gravitational binding energy of the star-planet pair is 
negative. Any object (e.g., a star, a binary system, a plane¬ 
tary system, or a flee-floating planet) is considered as having 
escaped from the star cluster when each of the three follow¬ 
ing conditions are satisfied: (i) its velocity is larger than the 
star cluster’s escape velocity at the location of the object; (ii) 
the velocity vector points away from the star cluster centre, 
and (hi) the object is located at a distance r > 2Rt from the 
cluster centre, where Rt is the cluster’s tidal radius (Eq. [^. 
All other objects are considered as gravitationally bound to 
the star cluster. The four categories and their abbreviations 
(BPS, BFP, EPS, and EFP) are listed in Table None of 
the planets in our models experiences a physical collision in 
our simulations, and therefore the total number of planets 
is conserved at any time: 

-NbPS + -NbPS + NbFP + -NeFP = Nplanet . (3) 


(e.g., Binney & Tremaine 19871. Here we have adopted a 
Galactocentric distance Dq « 8 kpc and a Milky Way-like 
galaxy with an enclosed mass of Mg = 5.8 x The 

initial tidal radius is Rt « 5.5 pc ~ 7rhm for the star clusters 
modelled in our study. 

All models are evolved for t = 50 Myr, the time beyond 
which planetary systems experience little further evolution. 
For a set of star clusters with an identical number stars, 
average stellar mass, and virial radius, the N-body time is 
proportional to the physical time scale. It should be noted. 


We define the corresponding fractions /bps, /bps, /bfp, 
and /bfp as the ratios between the number of planets in 
each category, relative to Apianet • Since each star is initially 
assigned a planetary companion and is initially bound to 
the star cluster, the initial fractions are /bps = 100% and 
/bps = /bfp = /bfp = 0%. After complete dissolution of a 
star cluster, /bps = /bfp = 0% and /bps + /bfp = 100%. 

A BPS can escape from the star cluster intact as an 
EPS, or it can be ionised in the star cluster after a close 
encounter with another star and become a BFP. These 


































4 Zheng, Kouwenhoven & Wang 



Figure 1. Three realisations of the star clusters at t = 0. The different panels show the position of the stars in the cluster for fractal 
dimensions (from left to right) D = 1.6, 2.3 and 3.0, respectively. 


BFPs gradually escape from the star cluster and become 
EFPs (Wanget al. [201^1. Capture of single stars into bina¬ 


ries (e.g., Goodman fc Hut||1993 Kouwenhoven et al.|2010 


Moeckel & Bate 2010 Moeckel & Clarke 20111, re-capture of 


BFPs (e.g., Parker fc Quanz||2012 Perets & Kouwenhoven 

2012) by stars, and exchange of planets between stars (e.g., 
Jilkova et al.|20l5 ) in the star cluster occur, but these pro¬ 
cesses are rare. Given enough time, the vast majority of the 
planets therefore follows one of these dynamical tracks: 


BPS EPS (track 1) 

BPS BFP EFP (track 2) 


( 4 ) 


As a result, dN-sps/dt ^ 0, dNpps/dt ^ 0, and dNppp/dt ^ 
0 at any time. In the thousands of simulations we carry out, 
dynamical binary star systems form, and several of these 
may host planets in (primordial) S-type configurations or 
(captured) P-type configurations. In our analysis, we con¬ 
sider these objects as BPSs or EPSs, depending on whether 
or not they are a cluster member. Although these are in¬ 
teresting processes that warrant further study, they barely 
affect the analysis of the planetary populations as a whole, 
and we will therefore not discuss these in this paper. 


3 RESULTS 

In this section we describe how the planetary population 
evolves different environments and how these results depend 
on the orbital parameters of the planets. The global evolu¬ 
tion of the star clusters is described in Section [3.11 In Sec¬ 
tion |3]^ and Section [3.3| we describe how the planetary pop¬ 
ulation evolves over time, and how this evolution depends on 
the initial semi-major axis of the planets. Subsequently, we 
study the dependence on initial amount of substructure in 
the cluster (Section |3.4[ ), and on the mass of the star clusters 
(Section |3.5[ ). 


3.1 Star cluster evolution 


The dynamical evolution of substructured star clusters has 
been studied extensively (see, e.g., Gartwright fc Whitworth] 


|2004[ [Goodwin fc Whitworth] |2004[ [Allison et al.||2009[ 



0 10 20 30 40 50 0 10 20 30 40 50 

Time (Myr) Time (Myr) 


Figure 2. The evolution of the half-mass radius r^m (^e/i) and 
the central mass density p {right) for D = 1.6, 2.3, 3.0 (averaged 
over the ensemble of realisations). The different curves in each 
panel indicate models with initial virial ratios of Q = 0.3 (black), 
Q = 0.5 (red) and Q = 0.7 (blue). 


Parker & Quanz 2012[ and references therein). As in pre¬ 
vious studies, we find that the initial spatial distribution of 
stars in each cluster (see, e.g.. Fig. changes drastically 
on a short time scale. These changes typically occur within 
several dynamical timescales, tdyn, which is usually defined 
as 


tdy 


2 X 10^ 




V IO6M0 


-1/2 


\ 3/2 

^hm \ 

I— 

1 pc / 


( 5 ) 


(e.g., Spitzer|[l987 ). Although our models do not satisfy all 
necessary conditions for applying this equation, particularly 
in the highly-substructured clusters, it can be used to es¬ 
timate a dynamical time scale of tdyn ~ 0.8 Myr in our 




















































Planetary systems in young star clusters 5 


reference model. After substructure has been removed and 
virial equilibrium achieved, the clusters slowly evolve due to 
two-body relaxation, stellar evolution, and the effect of the 
Galactic tidal field. Two-body encounters typically happen 
within a relaxation time. 


trlx — 


0 . 206 jVsr;^^^ 

PUMlPnA ’ 


( 6 ) 


where G is the gravitational constant, Na is the number of 
stars, and In A « In Ns is the Coulomb logarithm ([Binney "IF 


Tremaine|1987|| Heggie fc Hut|2d03l. Although the modelled 


star clusters are initially substructured, Eq. provides a 
reasonable first-order estimate for the initial relaxation time, 
which is roughly trix ~ 18 Myr. It must be noted, however, 
that substructured star clusters tend to evolve faster (see 
Section[2| and show earlier signs of mass segregation | Allison 
^[2C 

ig-ll 


et al. 12009 12010 \ than initially smooth star clusters. 


rhn 


In Fig. we show the evolution of the half-mass radius 
and the central mass density p of each star cluster. Fol¬ 


lowing Parker & Quanz (20121 we define the central density 


as the average stellar density within the half-mass radius 

^hm 5 i.6.5 


3Mci 


(7) 


Although this definition of the density allows us to describe 
the overall evolution of the star clusters, it is still an average, 
and therefore tends to overestimate or underestimate of the 
local stellar densities, particularly in highly-substructured 
star clusters (see, e.g., Parker fc Dale|[2013 [Parker et al. 


2014 Parker|2014 1. This results in a higher destruction rate 


as compared to estimates from Eq. particularly in the 
central region; we will address this issue in Section [3.4| 

At t = 0 Myr all star clusters have rhm = 1 pc, as 
specified by the initial conditions (note that our clusters 
are typically a factor 1.26 larger than those of [Parker &| 
[Q uanz | [M2| ). The time at which the star clusters experi¬ 
ence a gradual transition from substructured systems into 
a smooth, slowly-evolving systems, is roughly proportional 
to Q, and depends only mildly on D. Apart from the cases 
with Q = 0.7, the clusters initially contract, after which 
rhm increases with time. All clusters with Q — 0.3 reach 
a central density higher than 700 Mqpc~^, which corre¬ 
sponds to an average stellar separation of 10^ AU. The clus¬ 
ters Q = 0.5 only reach a modestly high central density of 
p ~ 200 starspc“^. The clusters with Q — 0.7 initially ex¬ 
pand, until a large number of stars reaches reaches the tidal 
radius and escape. Due to the large number of escapers for 
the Q = 0.7 models, cluster membership is decreased, and 
the bound stars remain in a cluster with a smaller half-mass 
radius and larger central density. For times t > trix, all clus¬ 
ters expand and decrease in central density. Those that are 
initially supervirial obtain a larger rhm and a lower p than 
star clusters starting out with Q = 0.5, while the star clus¬ 
ters with Q = 0.3 generally have intermediate values for rhm 
and p. 

The close encounter probability is roughly proportional 
to p, which changes by almost two orders of magnitude dur¬ 
ing the first 50 Myr, although local density variations due 
to initial substructure may differ substantially. The vast ma¬ 
jority of close encounters and planetary disruptions there¬ 
fore occur during the high-density phase during the first 



Time (Myr) 


Figure 3. The evolution of A^bpSi ATbePi and A^gpp in 

the star clusters (averaged of over the ensemble of realisations). 
All planets initially have a = fOO AU. Results are shown for star 
clusters with initial substructure parameters D = 1.6 (black), 
D = 2.3 (blue) and D = 3.0 (green) and initial virial ratios 
Q = 0.3 (solid curves), Q = 0.5 (dotted curves) and Q = 0.7 
(dashed curves). 


~ 20 Myr. Although initial substructure and the initial virial 
state mostly affect the evolution of the cluster at early times, 
they play an important role in determining the distribution 
of the planets over the four dynamical states at later times. 


3.2 Evolution of the planetary population 

3.2.1 Properties of the planet population 

The evolution of A^bps, A^bfp, A^eps and Abfp is shown 
in Fig. [^ for clusters with different D and Q. The plane¬ 
tary population experiences most evolution during the first 
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Figure 4. The velocity distribution of BPS (purple), BFP (red), EPS (blue) and EFP (black), at t = 50 Myr, for the combined sets of 
twenty realisations. 


~ 5 Myr. The most pronounced changes are seen for highly- 
substructnred clnsters and those out of virial equilibrium, 
resulting in smaller A^bps and larger AIbps, Nbfp and AIefp 
at any time. During the first ~ 5 Myr, A^bps decreases and 
AIbfp increases, indicating that many of the ionised planets 
initially remain part of the cluster. A substantial fraction of 
these BFPs have speeds up to a few kms“^ above the local 
escape velocity (see below). As planets are marked as esca¬ 
pers beyond two tidal radii, it takes St ~ 2Rt/v ~ 2 —Q Myr 
to escape. At later times, after the star clusters have ob¬ 
tained a state of quasi-equilibrium, Abps and Abfp both 


slowly decrease with time due to gradual disruption of BPSs 
and gradual escape of BPSs and BPFs. The nearly flat dis¬ 
tribution A^BFp(t) indicates that the production and escapes 
rate of BFPs are roughly equal for t > 10 Myr. As dynam¬ 


ical capture of free-floating planets is rare (e.g., Parker & 
|Quanz|[2012[ [Perets Kouwenhoven||2012[ ), the increase in 
Abps is almost entirely due to escape of existing planetary 
systems (Eq. |^. As virialised star clusters tend to lose stars 
at a roughly constant rate (e.g., [Heggie fc Hut||2003| ), and 
since none of the planetary systems is disrupted after escape, 
A^bps increases more or less linearly with time. 
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As expected, the dynamical evolution is strongest for 
the models with (D = 1.6, Q = 0.3) and weakest for those 
with {D = 3.0, Q = 0.7). At t = 50 Myr, the {D = 1.6, 
Q = 0.3) clusters have lost 35% of their BPSs: 4% have 
escaped the cluster intact as EPS, 12% escaped from the 
cluster as EFPs, and 19% remain part of the cluster as 
BFPs. On the other hand, at t = 50 Myr, the [D = 3.0, 
Q = 0.7) clusters only lose 11% of their BPS: 7% become 
EPSs, 2% remain in the cluster as BFPs, and 2% escape as 
EFPs. Clusters with other initial conditions show similar or 
intermediate behaviour. Variations in D and Q result in a 
scaling down of A^bps, and in a scaling down of the equilib¬ 
rium value (beyond t > 10 Myr) for Nbfp, and also results in 
an upward scaling of Neps and Vefp- An EFP is generated 
if and only if a BFP escapes. As free-floating planets tend 
to escape from star clusters substantially earlier than stellar 
members (e.g., Parker fc Quanz]|2012 Wang et al. 2015aI, 
Abfp increases fast at early times, after which it grows at 
a steady rate when the BFP equilibrium is established. As 
all clusters are smooth and virialised beyond t « 10 Myr, 
the differences at later times can be explained as being the 
result of the initial dynamical evolution of the star clusters 
under different values of D and Q, more specifically, as a 
result of the variations in stellar density during the early 
phase of evolution (see, e.g., Parker et al.||20lTK 


Fig .[^ shows the velocity distributions at 50 Myr. Those 
of the BPSs and BFPs are more or less identical, although 
BFPs tend to roam around in the outskirts of the cluster and 
these therefore have slightly lower velocities. The escapers 
(EPS and EFP) have higher velocities, usually 1 — 3 kms“^ 
above the escape velocity, consistent with the findings of 


Parker & Quanz (20121. Some EFPs reach velocities as high 
as ~ 30 kms“^ and escape as ” runaway planets” following a 
strong dynamical interaction. There is no strong correlation 
between the velocity distribution of the EPSs and the initial 
conditions. 


Non-disruptive dynamical encounters tend to alter the 
orbital elements of BFSs, as shown in |Spurzem et alT] ( |2009[ ) 
and Parker & Quanz (20121. In Table we summarise the 
semi-major axis a and eccentricity e distributions for our 
clusters at t = 50 Myr. Weak encounters tend to result in 
more angular momentum exchange (changes in e and i) than 
energy exchange (changes in a), resulting in the familiar 
fountain-diagrams (e.g., figure 5 in [Parker fc Quanz||2012] ). 
Eccentricity growth is largest for highly substructured, sub- 
virial clusters, where dynamical interactions are frequent. 
Semi-major axes remain mostly unchanged, although the 
fraction of softened BPSs (larger a) is again largest for clus¬ 
ters with D — 1.6 and Q = 0.3, as close encounters in these 
clusters are most frequent. On average, 70—90% of the BPSs 
that remain at t = 50 Myr have slightly hardened (smaller 
a), and this fraction is largest for the least violent clusters. It 
should be noted, however, that most of the hardened orbits 
experience a minimal decrease in semi-major axis, while the 
softened orbits can become substantially wider (see, e.g., hg- 
ure 5 in |Parker fc Quanz|2012| and note that the horizontal 
axis is logarithmic). 


3.2.2 The fraction of planet-hosting stars in star clusters 
and in the field 


In the previous sections we described the dynamical prop¬ 
erties of the entire set of planets. Dynamically, this makes 
sense, as we can follow all planets in our simulations. Obser- 
vationally, it is more useful to consider the fraction of planet¬ 
hosting stars in the cluster and the held, and the number 
of free-hoating planets in the cluster and in the held, with 
respect to the number of stars in the same environment. At 
a given time, the fraction of planet-hosting stars /p,b in the 
star cluster and /p,e in the held are 


7p,b 


Abps 

Sb Mbps 


and /p,e 


Veps 

Se NbpS 


( 8 ) 


respectively. Here, Sb and Se represent the number of stars 
without a planetary companion in the star cluster and in 
the held, respectively. Similarly, the ratio between the num¬ 
ber of free-hoating planets and stars among the bound and 
unbound objects is 


/ff.b = 


Vbfp 
Sb -|- Nbps 


and /ff,e = 


Aefp 

Se -f VePS 


(9) 


Since all planets initially orbit a cluster member, the ini¬ 
tial values are /p.b = 1, /ff,b = 0, while /p,B and fs^e are 
undehned. 

Fig. shows /p,b, /p.e, /ff.b, and /ff,e as a function of 
time for clusters with different initial Q and D. A compar¬ 
ison between Figures and shows a strong correlation 
between /p^b and Vbps, which is not surprising since the 
number of stars without a planet in the star cluster {Sb in 
Eq. 1 ^ is generally much smaller than Nbps. During the ini¬ 
tial virialisation process the fraction /p^b rapidly decreases, 
after which it obtains its equilibrium at a value that depends 
on D and Q. The ratio /p,e for the escaped stars is unde¬ 
hned during the hrst ~ 5 Myr, as no star has had the time to 
escape from the cluster. Subsequently, /p,e grows with time 
(after a short period strong huctuations due to low-number 
statistics), because BFPs tend to escape at earlier times than 
EPS (see Fig.j^ as on average they obtain higher velocities 
after experiencing close encounters with other cluster mem¬ 
bers. The fractions /p,b and /p_b show a similar dependence 
on D and Q. 

As Sb is small compared to Abps at any time, the evo¬ 
lution of /fi,b also correlates strongly with /bps (Fig. [^. 
The ratio /ff,e among the escaping objects is initially unde¬ 
hned, and as initially only free-hoating planets escape from 
the cluster, the value is initially above unity. At later times, 
stars also escape and the values drop below unity. For the 
models with Q = 0.7 the ratio /ff,e drops to low values (more 
stars than EFPs), as most planetary systems escape intact 
(as EPSs), while for highly substructured and subvirial star 
clusters, /ff,e remains above unity (more EFPs than stars) 
at t = 50 Myr. 


3.3 Dependence on initial semi-major axis a 

The survival chances of planetary systems depend on both 
their environment and on their initial binding energy, Eb oc 
a~^. We study the latter dependence by modelling the evo¬ 
lution planetary systems over a large range of semi-major 
axes, 1 AU a ^ 10 000 AU, and study their evolution for 
an ensemble of eight hundred star clusters. Fig.j^shows the 
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Table 3. The semi-major axis a and eccentricity e of the remaining planetary systems at t = 50 Myr, for star clusters with different initial 
substructure parameters D and virial ratios Q. The ’’ps” symbol indicates changes of less than one part per million. The percentages 
represent the averages of twenty realisations. The range of values for the ensemble of simulations is indicated between the brackets. 


D 

Q 

a 

100 AU, e Ri 0 

a 

100 AU, e > 0 

a > 100 AU, all e 

a < 100 AU, all e 

all 

a, e ^ 0.1 




per cent 


per cent 


per cent 


per cent 


per cent 

1.6 

0.3 

4.4 

(0.0 - 28.6) 

4.9 

(1.3 - 21.7) 

17.7 

(12.0 - 22.7) 

72.6 

(30.1 - 84.7) 

24.7 

(17.6 - 32.4) 

2.3 

0.3 

0.6 

(0.0 - 1.5) 

2.5 

(0.8- 5.7) 

15.6 

(13.4 - 20.3) 

81.2 

(77.9 - 85.8) 

19.4 

(16.4 - 24.4) 

3.0 

0.3 

0.7 

(0.0 - 1.7) 

2.2 

(0.9-4.0) 

14.6 

(11.0 - 18.1) 

82.5 

(77.5 - 87.5) 

17.6 

(12.6 - 22.5) 

1.6 

0.5 

0.7 

(0.0 - 1.6) 

1.8 

(0.6- 5.0) 

12.3 

(9.3 - 17.8) 

85.1 

(80.4 - 88.8) 

17.3 

(12.0 - 22.6) 

2.3 

0.5 

0.9 

(0.2 - 1.6) 

1.2 

(0.2 - 3.3) 

10.3 

(7.3 - 13.3) 

87.6 

A-6 - 90.6) 

12.9 

(9.4 - 16.8) 

3.0 

0.5 

0.7 

(0.2 - 1.4) 

1.2 

(0.0- 2.0) 

10.4 

(8.5 - 16.4) 

87.7 

(81.4-90.7) 

11.4 

(7.3 - 14.6) 

1.6 

0.7 

1.1 

(0.0 - 7.2) 

1.5 

(0.3-4.9) 

10.2 

(6.0 - 14.2) 

87.2 

(73.6 - 93.7) 

13.6 

(9.4 - 19.1) 

2.3 

0.7 

0.8 

(0.2 - 1.6) 

0.8 

(0.0- 2.0) 

9.0 

(6.7- 12.4) 

89.4 

(84.2 - 91.8) 

9.1 

(4.7-11.9) 

3.0 

0.7 

1.1 

(0.2 - 2.4) 

1.2 

(0.2 - 2.0) 

8.6 

(5.9 - 11.2) 

89.1 

A-5 - 91.8) 

7.4 

(5.8 - 9.5) 
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Figure 5. The fraction of stars with planets /p^b star cluster {top-left)^ the fraction of stars with planets /p,e among the escaping 

stars (bottom-left), the ratio between free-floating planets and stars in the cluster /ff,b (top-right), and the ratio between free-floating 
planets and stars among the escaped objects /ff g (bottom-right) as a function of time (Eqs. [^and[^. Results are shown for star clusters 
with initial substructures D = 1.6 (black), D = 2.3 (blue) and D = 3.0 (green) and initial viral parameters Q = 0.3 (solid curves), 
Q = 0.5 (dotted curves) and Q = 0.7 (dashed curves). 
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Figure 6. The frequency of BPS (black), BFP (red), EPS (blue), and EFP (green) at f = 50 Myr, as a function of initial semi-major 
axis. The black and red curves are well fitted by the distribution f{a) = /o(l H- [n/ao]*^)”^, and the fitted values are listed in Table[^ 


Table 4. The fitted values for the distribution /o(l -l- [a/ao]^) ^ 
for bound planetary systems and bound free-floating planets in 
Fig. at t = 50 Myr. 


D 

Q 

/o 

BPS 

ao 

c 

fo 

BFP 

ao 

c 

1.6 

0.3 

83.3 

226.0 

1.77 

91.3 

235.9 

-1.69 

2.3 

0.3 

90.1 

317.2 

1.94 

92.7 

334.5 

-1.95 

3.0 

0.3 

92.0 

348.5 

2.02 

92.9 

370.0 

-2.07 

1.6 

0.5 

86.8 

270.7 

1.81 

93.3 

283.5 

-1.73 

2.3 

0.5 

93.0 

371.5 

1.97 

94.2 

393.9 

-2.03 

3.0 

0.5 

93.8 

405.7 

2.02 

94.3 

427.2 

-2.10 

1.6 

0.7 

84.8 

315.3 

1.85 

89.3 

322.9 

-1.74 

2.3 

0.7 

90.8 

426.1 

2.03 

90.6 

440.4 

-2.13 

3.0 

0.7 

92.0 

454.0 

2.02 

90.5 

470.5 

-2.17 


dynamical fate of planets with different semi-major axes at 
t = 50 Myr for the different clusters. 

Most planetary systems with a < 100 AU survive as star 
cluster members, although a small fraction escape as EPSs. 
Almost all systems with a > 2000 AU are destroyed and the 
majority of these remain as BFPs in the star cluster, while 


some escape as EFPs. Survival chances drop sharply from 
a ~ 100 AU to a ~ 2000 AU, and as a result, /bfp increases 
and /bps decreases rapidly with increasing a. As discussed 
earlier (see Fig. |^, most of these changes occur during the 
hrst ~ 5 Myr. The fraction of planetary systems that escape 
the cluster intact (EPS; blue symbols) is typically 5 — 15% 
when a < 200 AU. The number of EFPs slowly increases 
with increasing a, as more systems are disrupted. The num¬ 
ber of EFPs also drastically increase around a « 200 AU. 

In order to quantify our results, we ht /bps (a) and 
/bfp (a) to distinguish between the effects of different ini¬ 
tial conditions in Fig. All curves are well described with 
a function 


/(a) 


/o 




( 10 ) 


where fo,ao and c are constants. For a <C oq, we ob¬ 
tain /bps (a) + /bps (a) « 1 and /bfp (a) = /bps (a) « 0, 
which refers to the initial conditions. For ao we obtain 
/BPs(a)+/BPs(a) ~ 0, indicating that almost all systems are 
destroyed. The fitted quantities /o, no and c of the BPS and 
BFP populations are listed in Table for models with dif¬ 
ferent initial D and Q. For comparison, the fits for /bps (a) 
and /bfp (a) are also combined in Fig. 

Strong transitions occur in the region around a = oq. 
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Figure 7. Fitted curves describing the frequency /bps 

/bfp (bottom) at t = 50 Myr, as a function of semi-major axis 

(cf. Fig.|^. 


which is the semi-major axis boundary that separates the 
domains in which planetary systems remain mostly intact or 
are mostly destroyed. The quantity |c| indicates how broad 
this region transition is in terms of semi-major axis. Note 
that a smaller value of |c| corresponds to a broader transition 
region. Note the value of no is unrelated to the hard-soft 
boundary that is often used to describe the dynamics of 
binary systems in star clusters; all planets are weakly bound 
to their host star and are therefore by definition ’’soft”. The 
destruction of planetary systems is determined by the close 
encounter rate, rather than by the encounter energy. 

For the BPSs, the quantity /o « 83 — 94% indicates the 
fraction of stars remaining in the star cluster. Since almost 
all planetary systems remain intact when a ao, f{a) « fo 
for the BPS and /(a) — 1 — fo for the EPS in this regime. 
In other words, 6 — 17% of the planetary systems escape 
the cluster within 50 Myr. As described earlier the number 
of planetary systems that remain part of the star cluster 
(/o) increases with increasing degree of homogeneity, and is 
largest for clusters starting out in virial equilibrium. In the 
limit a flo, virtually all planetary systems are destroyed, 
and free-floating planets either remain in the cluster as BFPs 
or escape as EFPs. In this limit, /(a) = fo ~ 91 — 94% 
for the BFPs, and /(a) = 1 — /o = 6 — 9% for the EPSs. 
This demonstrates that even when planetary systems are 



Figure 8. The dependence of fBPs{D) on the initial substructure 
parameter D. Shaded regions bracket, from dark to light, /bps W) 
at t = 10, 20, 30, 40, and 50 Myr. The dots represent the median 
values (and corresponding standard deviations) at t = 50 Myr, 
and the solid curve represents a fit to these data (Eq. [T^. 


destroyed at early times, most free-floating planets remain 
a member of their host star cluster during the first ~ 50 Myr. 
The critical semi-major axis oo is in the range 200 — 500 AU, 
with the exact value depending on the initial conditions. 
BPSs in more substructured star clusters experience more 
close encounters, and therefore have a smaller ag. Planetary 
systems have the largest chances of being retained (i.e., no 
disruption or escape) in the least violent clusters (Q = 0.7 
and D = 3.0), and these clusters thus have the largest ag. 

The width in the transition region is described by the 
parameter |c|, where the sign of c determines whether the 
distribution increases or decreases with a. Part of this vari¬ 
ation is a result of the global density variations in the star 
clusters, and part can be attributed to the variation in bind¬ 
ing energies (host stars with different stellar masses). For 
the fits of both the BPS and the BFP, |c| increases with 
increasing D, indicating that the transition range (around 
ag) for destruction of planetary systems is broader. This 
is understandable, since an increased amount of initial sub¬ 
structure results in a larger variation in local stellar densities 
and therefore also encounter frequencies. For the BPSs, |c| 
depends only mildly on Q, while for the BFPs, |c| increases 
with increasing D, for the same reason as above. 

Bonnell et al. | pooI| obtained comparable results using 


analytical estimates (see their figures 4 and 5). They study 
the evolution of planetary systems in different environments 
over much longer time scales (up to 10 Gyr), and focus on 
planetary semi-major axes in the range 0.1 — 100 AU. For 
the cluster densities considered in our work, only the widest 
planetary orbits in this separation range are affected at t = 
50 Myr (Fig. [^, which is indeed what [Bonnell et ^ ( |2001[ ) 
have found. Further V-body simulations spanning a much 
larger range in cluster density and a much longer integration 
time are required to validate the other analytical results 
presented in Bonnell et al. (20011. 
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3.4 Dependence on the substructure parameter D 

The encounter rate T{D) of stars in a star cluster with sub¬ 
structure parameter D can, to first order, be approximated 
with r(D) = n{D)'Kp^^^v{D), where, Pmax is a maximum 
impact parameter that is considered as encounter, n is the lo¬ 
cal stellar density and v{D) is the relative velocity at inhnity 
between two bodies. For highly substructured star clusters, 
n{D) and v{D) vary strongly with location due the presence 
of high-density pockets of stars, a phenomenon that was dis¬ 
cussed earlier by, e. g., [Parker fc Dale| ( |2013[ ); [Parker et al.| 
(20141 and Parker 120141. In those cases, a larger number 


of planetary systems is disrupted, resulting in a larger /bps 
and a smaller /bfp. 

To quantify the correlation between D and /bps, we run 
simulations with thirty different initial substructure param¬ 
eters in the range D £ [1.6 — 3.0]. All cluster are initially 
in virial equilibrium {Q = 0.5), and the remaining initial 
conditions are identical to those of our reference model (Ta¬ 
ble [^. In order to reduce statistical fluctuations, we run 
twenty realisations for each value of D. The resulting values 
fBPs(D,t) are shown in Fig. for t = 10, 20, 30, 40, and 
50 Myr. As can also be seen in Fig.[^ /bps decreases roughly 
linearly with time, after the initial phase of relaxation when 
substructure is removed and virial equilibrium is restored. 
For the clusters in our sample, the fraction /bps(D, t) is well 
described by 


fBFP {D, t) ^ A{t) + B{D) 


( 11 ) 


where the (nearly linear) evolution over time can be ex¬ 
pressed as 


A{t) « 1 - (2.06 X 10"^)t-b (3.28 x W~^)P 


( 12 ) 


where t is in units of Myr. The dependence on the initial 
substructure parameter as 

B(D) « (0.78-(13) 

The approximation in Eq. El has of one term that depends 
only on t and another that only depends on D. This means 
that changes over time, dfBPs/dt = dA{t)/dt, are mostly in¬ 
dependent of the initial amount of substructure, which can 
be seen in Fig. It also means that a different amount of 
initial substructure does not affect its evolution over time, 
dfsps/dD = dB{D)/dD, but merely establishes its normali¬ 
sation at later times. At any moment in time beyond several 
initial relaxation times, /bps depends weakly on the initial 
amount of substructure for D > 2, when d/pps/dD ~ 0, 
which explains the similarities in Fig. Note that, as ex¬ 
pected, /BPs(D,t = 0) « 1 corresponds to our initial condi¬ 
tions for D > 2. For D < 2, /bps(D, t = 0) > 1, which shows 
that our simple approximation is not valid during the earli¬ 
est phases of evolution of star clusters with a high amount 
of initial substructure. 


3.5 Dependence on initial mass and density 

In the previous sections we discussed the evolution of star 
clusters that all consisted of an identical number of stars 
(As = 1000) and had an identical size (Phm = 1 pc). In 
order to broaden the applicability of our results, we discuss 
in this section how the evolution of the planetary population 



Number of stars 


Figure 9. The frequencies /bps, /bps, /bfp and /efp at t = 
50 Myr, as a function of the number of stars in the cluster As. 
Results are shown for two clusters with central stellar densities 
(Eq.[^ Pi = 120 pc“^ and p 2 = 280 pc“^. 


depends on the initial mass and central density of the star 
clusters. 

We carry out simulations with As = 25, 50, 100, 250, 
500, and 1000 stars. In order to allow a fair comparison be¬ 
tween the clusters of different mass, one would ideally keep 
all other star cluster parameters (such as its half-mass ra¬ 
dius rhm, its density p, and its half-mass relaxation time 
frix) constant. As all these parameters depend on As and 
on each other, only combinations of parameters can be kept 
constant (see jKouwenhoven et al.|2014 for an extensive dis¬ 
cussion on how different combinations of these parameters 
affect star cluster evolution). As we are mostly interested 
in dynamics, we compare star clusters in which the initial 
stellar density is constant for all values of As. We model two 
sets of star clusters: one set in which all star clusters have 
central stellar densities (Eq. of pi = 120 pc“® (which in¬ 
cludes our reference model) and another set in which all clus¬ 
ters have p 2 = 280 pc“®. Consequently, the corresponding 
half-mass radii can be expressed as rhmi = O.lAs^^® pc and 
nim 2 = 0.075As^^® pc. All star clusters initially have D = 3.0 
and Q = 0.5, and all stars initially host a planet in a circular 
orbit at a = 100 AU. The fraction of planets in the different 
dynamical categories at t = 50 Myr are shown in Eig.[^ and 
results are obtained by averaging ten realisations for each 
model. 

The (initial) dynamical time scales as fdyn oc p~^^^ 
(Eq. [^, and is therefore identical for star clusters over a 
large mass range, with a given initial stellar density, corre¬ 
sponding to tdyn ~ 0.82 Myr for pi and 0.53 Myr for p 2 . The 
relaxation time, on the other hand, scales approximately as 
trix oc (Eq.|^, and ranges between 1 Myr to 18 Myr 

for the least and most massive clusters in our sample, respec¬ 
tively. The dynamical age of the star clusters at t = 50 Myr 
is thus proportional to Aix oc tdyn/Aa and is independent 
of p. Finally, the tidal radius (Eq. of a cluster scales as 
Rt oc A^/®, and is also independent of p. 

To understand dependence on the initial density in 
Fig.ii it is convenient to first consider the evolution of the 
stellar population alone, and ignore the presence of plan- 
































12 Zheng, Kouwenhoven & Wang 


ets, which have a negligible effect on the stellar dynamics. 
Clusters with a given are all evolved for an identical 
number of relaxation times, despite the different stellar den¬ 
sities. The only relevant difference between the clusters is 
the tidal radius, more specifically the ratio rhm/Rt cx p~^. 
Clusters with a larger initial density are more compact, and 
it therefore takes longer for BPSs and single stars to escape 
while the clusters gradually expand to fill their tidal radii. 
Consequently, clusters with larger p have a larger number of 
bound stars, after an identical number of dynamical times 
have elapsed. Larger densities, however, do result in the de¬ 
struction of more planetary systems, and therefore exhibit 
a larger fraction of BFPs and EFPs at 50 Myr. 

Fig.|^ shows that /bps increases with increasing star 
cluster mass. In general, /bps on both the stellar density, 
and the binding energy of the star cluster. The former is 
related to the strength and frequency of close encounters, so 
for the two different star density cases, it is reasonable to 
see that when a cluster has a higher initial density, it has 
smaller /bps at t = 50 Myr. Due to the low cluster masses, 
many of the liberated planets are ejected from their host 
stars with velocities above the cluster’s escape velocity, and 
therefore /bfp is close to zero and /efp is positive (apart 
from the higher-mass clusters that can retain free-floating 
planets). Planetary systems escaping the cluster intact are 
most common among the lowest-mass clusters that experi¬ 
ence near-dissolution around t = 50 Myr. 


4 CONCLUSIONS AND DISCUSSION 


Although many surveys have focused on detecting planets 
in open and globular clusters, relatively few planets orbit¬ 
ing stars and free-floating planets have been discovered in 
these environments. This can partially be attributed to ob¬ 
servational difficulties, but also to disruption following close 
encounters with neighbouring stars and subsequent escape 
of free-floating planets from star clusters. Inspired by the 
earlier study of Parker & Quanz (20121 we carry out N- 
body simulations of evolving star clusters. We employ the 
NB0DY6 package to model open star cluster environments 
with Ns = 1000 stars, in which we assign each star a Jupiter- 
mass companion. All planets initially have circular and ran¬ 
domly oriented orbits. We evolve all clusters for 50 Myr and 
study how the evolution of the planetary population depends 
on the initial conditions. We study star clusters with differ¬ 
ent initial amounts of substructure {D = 1.6 —3.0), different 
initial virial ratios {Q — 0.3 —0.7) and initial planetary semi¬ 
major axes in the range 1 AUsi a ^ 10 000 AU. To broaden 
the applicability of our results, we also carry out simula¬ 
tions of star clusters with different total masses and stellar 
densities. Our main results can be summarised as follows: 


(i) The initial values of D and Q mostly affect the evo¬ 
lution of planetary systems during the first ~ 5 Myr, which 
is the time at which the initial substructure is removed and 
virial equilibrium is achieved. Subsequently, the stellar and 
planetary populations evolve gradually, more or less inde¬ 
pendent of the initial conditions. Although the initial choices 
for D and Q take effect at early times, their influence on the 
fate of planetary systems can be substantial. The amount of 
disruption of planetary systems and escape of free-floating 
planets is larger when a star cluster is initialised with a 


larger amount of substructure or is further away from virial 
equilibrium. 

(ii) Although both play a role in the disruption rate of 
planetary systems and the escape rate of planetary systems 
from the star cluster, the initial amount of substructure, D, 
mostly affects the former, while the initial virial state, Q, 
has most influence on the latter. 

(iii) In addition to environmental factors, the disruption 
rate of planetary systems is strongly correlated with initial 
semi-major axis a. The fractions of bound planetary systems 
/BPs(fl) and bound free-floating planets /BFp(n) are well- 
described by the functional form /(a) = /o(l -I- (a/oo)'’)”^, 
where /o, oo and c are constants. In the case of /bps, /o rep¬ 
resents the fraction of stars that remains bound the cluster, 
ao is the stability limit for disruption of planetary systems, 
and c measures the width of the transition region, which is 
negative for /bps (a) and positive for /bfp(o). 

(iv) A higher degree of initial substructure in the star 
clusters results in a higher disruption rate of planetary sys¬ 
tems. The fraction of bound planetary systems over time is 
well approximated with /BPp(D,t) = A(t) + B{D), where 
A{t) is an almost linear function of time (Eq. 121, and B{D) 
is a function that increases with D and flattens off when the 


level of substructure approaches D = S (Eq. 131. 


(v) For clusters with a fixed initial density, the fraction 
of planetary systems present at t = 50 Myr increases with 
increasing cluster mass, while the fraction of escaping free- 
floating planets decreases. Crowded environments result in 
more frequent encounters, so more free-floating planets are 
generated in the high-density centre of the cluster, and fewer 
intact planetary systems are ejected from the star cluster. 


Our study has focused on the general dynamical evolu¬ 
tion of star-planet systems and free-floating planets in star 
clusters. Our models represent a simplification of reality, as 
we have exclusively modelled single-planet systems. Multi¬ 
planet systems are substantially more fragile than single¬ 
planet systems. Small perturbation of an outer planet can 
induce strong gravitational interactions with other planet (s) 
in the system, which may result in a reconfiguration of the 
system, in the ejection of one or more planets, or in a phys¬ 


ical star-planet or planet-planet collision (see, e.g., Hao et 
aI.|2013|[Shara et al.|20l4 |. Moreover, we have not included 


primordial binaries. Observations have indicated that stars 
of all masses and ages are often part of a binary or multiple 


stellar system (see, for example. 

Duquennoy & Mayor 

1991 

Kouwenhoven et al.||2005l |2007 

Connelley et al.| 2008a|b 

Raghavan et al.||2010 Bergfors et al.||20101 |Janson et al. 

2014 

Tokovinin|2014ajb and numerous others), and that a 


considerable fraction of the known exoplanets are part of a 
multi-planet system (e.g., Latham et aI.|20TT|. The presence 


of binary and multiple stellar systems increases encounter 
rates due to their larger collisional cross-section, and also 
extends the longevity of star clusters. Although dynamical 
binary systems form in our modelled star clusters, primor¬ 
dial binary systems are not included at this stage. 

The simplifications mentioned above warrant a further 
dynamical study in which stellar and planetary multiplicity 
is taken into account. It is technically very difficult (but not 


impossible; see 

Malmberg et al.||20111 |Davies et al.||20131 

|Cai et al.|2015b 

1 to model the evolution of multi-planet sys- 


terns in star clusters due to the enormous dynamical ranges 
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in size, time, and mass. However, with the recent upgrade 
of NBODY6++ (NBODY6++GPU; see |Wang et al.|2015b l 
and its integration into the AMUSE framework (Portegies 


Zwart et al. 2013[ Pelupessy et al.|2013[ Cai et al.|2015a I, di 
rect A^-body simulations of multi-planet systems in massive 
star clusters may be carried out in the near future. 
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